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"^J ■ Molecular dynamics algorithms are subject to some 
_— , amount of error dependent on the size of the time 
<0> step that is used. This error can be corrected by 
periodically updating the system with a Metropolis 
Q_priteria, where the integration step is treated as a 
selection probability for candidate state generation. 
Such a method, closely related to generalized hybrid 
Monte Carlo (GHMC), satisfies the balance condi- 
Cn tion by imposing a reversal of momenta upon can- 
didate rejection. In the present study, we deraon- 
strate that such momentum reversals can have a sig- 
Q-mificant impact on molecular kinetics and extend the 
Q_time required for system decorrelation, resulting in 
an order of magnitude increase in the integrated auto- 
g correlation times of molecular variables for the worst 
q cases. We present a simple method, referred to as 
^ reduced-flipping GHMC, that uses the information of 
Q the previous, current, and candidate states to reduce 
* ^5 * ne probability of momentum flipping following candi- 
>">date rejection while rigorously satisfying the balance rrru 
*£h condition. This method is a simple modification to neory 
(^^traditional, automatic-flipping, GHMC methods and 
significantly mitigates the impact of such algorithms 
,— ' on molecular kinetics and simulation mixing times. 

^ 1 Introduction 
in 

On 



and GHMC methods for Langevin dynamics with a low fric- 
tion constant result in a dramatic increase in the time required 
for system decorrelation. 

Recent results have demonstrated that an implementation 
of GHMC with no momentum flipping results in simulation 
data that closely reproduce the desired distribution, though 
the algorithm cannot be proven to satisfy the balance con- 
dition that ensures proper sampling As we demonstrate 
below, these results suggest a level of quantitative similarity 
for the average acceptance rates for two states with identical 
coordinates but opposite momenta. Though this equivalence 
is not generally true, we propose a method that uses a simu- 
lation's immediately available information-that of the previ- 
ous, current, and candidate states-to reduce the rate at which 
momenta are nipped. This method represents a simple mod- 
ification to the GHMC algorithm that is generally applicable 
to integration schemes and provides a significant benefit to 
simulation mixing times. 



2.1 Traditional GHMC 

Let x 6 f2 denote some state of our system observed with 
probability 

P(x) = iexp[-/3tf(x)] (1) 



* Molecular dynamics integration algorithms can be subject to 



Z = 
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some amount of error resulting from finite time step effects, 
£vq resulting in sampling that deviates from the desired station- 
i-H ary distribution. These effects can be corrected using hy- 
^ brid Monte Carlo (HMC) [l| or generalized hybrid Monte 
Carlo (GHMC) Q algorithms that implement a Metropolis 
acceptance criteria in which the dynamics paths are treated 
^ as candidate selection probabilities. If the momenta are not 
" randomized at the refreshment steps (as is done in the tradi- 
tional HMC algorithm) the detailed balance condition for this 
method imposes that the algorithm either reverse momenta 
upon rejection or that it reverse momenta upon acceptance 
(indicative of the GHMC algorithm) [1(3]. 

The momentum adjustments required by these algorithms 
elicit concern that their implementation will diminish sim- 
ulation mixing times or encumber interpretation of system 
dynamics [3|-|5[ . The results of the present study demonstrate 
that this is indeed the case: figure [1] shows that, even with 
a small timestep and high acceptance rate, use of the HMC 



[ dxexp[-/3i? (x) 
Ja 



(2) 



where j3 = is the inverse thermal energy. Also note that 
P (x) = P (x) , the tilde denotes momenta- reversal, and the 
Hamiltonian H is invariant under such a transformation. We 
will assume a general simulation framework, where a candi- 
date state y is generated from current state x with an associ- 
ated selection probability S (x, y) . It is assumed that y is gen- 
erated using an integration method that is irreducible (signi- 
fying that it can reach any state of the system from any other 
state), necessary to ensure ergodicity once Metropolization is 
included. For the overall transition probability T (x, y) , the 
detailed balance condition imposes that 



P(x)T(x,y) =P(y)T(y,x) 
The Metropolis criterion satisfies this condition 

5(y,x) ^(y) 

S(x,y) P(x) 



(x,y) = min 1, 



(3) 



(4) 
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There are a number of methods appropriate for generat- 
ing candidate state y in this framework, including Langevin 
integrators integration using the Andersen [ll| or 
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Figure 1: Comparison of the autocorrelation function for a 
polymer's radius of gyration obtained from simulations with no 
Metropolis test and with the GHMC and HMC algorithms. Sim- 
ulations were performed with a 10 fs time step and a friction con- 
stant of 7 = 0.01 ps _1 . See methods section for full simulation 
details. 



stochastic velocity rescaling thermostat [12j, and the tradi- 
tional GHMC algorithm (leap-frog integration coupled with 
a mixing angle that introduces random velocity collisions) Q 
. Simulations of this study generate candidate states with a 
single integration step using the velocity verlet with veloc- 
ity randomization (VVVR) integrator [lOj, equivalent to the 
Langevin integration scheme of Bussi and Parinello or 
the stochastic velocity rescaling thermostat 12|. It should 
be noted that the stochastic velocity rescaling thermostat 
appears to exhibit irreducible behavior (l2l |. but this has 
not been rigorously proven. Selection probabilities for the 
stochastic velocity rescaling thermostat are given in appendix 
[S]and can be obtained for the VVVR integrator using the ab- 
solute path action derived in reference |10| . 

We wish to demonstrate that, when coupled with either 
momentum reversal upon rejection or reversal upon accep- 
tance, equation [3] satisfies the balance condition P (x) = 
f n tiyP(y)5(y, x)P A (y,x). This can be shown following a 
previously outlined proof [l3l |. for which we denote the aver- 
age acceptance rate for transitions initiated from x: 



The balance condition is also satisfied if the algorithm imple- 
ments momentum flipping upon move acceptance rather than 
rejection 



14] 



r(x) 



J dyS(x,y)P A (x,y) 



(5) 



If we begin in state y which we assume to have been drawn 
from the appropriate distribution with probability P(y), the 
next value in our chain, x, is drawn from a probability P* (x) 
with two contributions: acceptance of the candidate x from 
initial state y, and rejection of some candidate initiated from 
y = x. Summing these contributions and using detailed bal- 
ance: 

P*(x) = ydyP(y)5(y,x)P A (y,x)+P(x)[l-r(x)] 

= J dyP(x)S(x, y)P A (x, y) + P(x) [1 - r(x)] 
= P(x)(r(x) + [l-r(x)]) = P(x) (6) 



Recent results have demonstrated that an implementation 
of GHMC with no such momentum flipping is capable of accu- 
rately reproducing certain statistical properties of molecular 
simulations This is a curious result, given that this al- 
gorithm cannot be shown to satisfy the balance condition. 
We note, however, that if some configuration x and the as- 
sociated x have identical acceptance rates, r(x) = r(x), then 
the balance condition is indeed satisfied without the need for 
momentum flipping: Following the proof given in equation [6l 
this can be shown assuming that r(x) = r(x) and summing 
the appropriate contributions: 

P* (x) = J dyP(y)S(y, x)P A (y, x) + P(x) [1 - r(x)] 

= P(x)(r(x) + [l-r(x)]) = P(x) (7) 

Though we generally expect that r(x) ^ r(x), the aforemen- 
tioned results Q suggest that r(x) and r(x) may have a level 
of quantitative similarity for some (perhaps easily met) set 
of conditions. Here, we develop a moveset that expands on 
this concept to use information from both the candidate con- 
figuration and the most recent configuration to reduce the 
probability of flipping momenta upon rejection while rigor- 
ously satisfying balance. 

2.2 Reduced-flipping GHMC 

For the traditional GHMC framework, Sohl-Dickstein Q has 
shown that, with proper accounting of the rates of inflow and 
outflow from x, the number of momentum flips upon candi- 
date rejection can be reduced while still ensuring satisfaction 
of the balance condition. This concept can be extended to 
algorithms that propagate coordinates and momenta in any 
fashion as long as the selection probability S is well defined. 

Our method, referred to as reduced-flipping GHMC, will 
select some candidate state (z) from the current state (y) 
according to the probability S (y, z). The acceptance proba- 
bility for this move, P A (y, z), obeys the standard Metropolis 
criterion given in equation 2] Upon rejection, we will define 
probabilities for flipping (y — > y) and self (y — > y) transitions 
that depend on both the candidate z and the previous state x. 
These transitions are associated with probabilities Pg(y|x, z) 
for self-transitions and Pp(y,y|x, z) for momentum flipping 
transitions. Thus, this algorithm proceeds as follows 

1. Select a candidate state z from current state y according 
to the probability S (y, z). 

2. Generate a uniform random number, £ £ [0, 1) and per- 
form the transition (y — > y2), where 



Y2 



if£<P A (y,z) 

ifP A (y,z) <£< (P A (y,z) 

else 



Ps(y|x,z)) 



where P A (y,z) +P s (y|x,z) +P F (y,y|x,z) = 1. 

The traditional GHMC method described in section 12.11 
henceforth referred to as the automatic-flipping algorithm, 
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always flips momenta upon candidate rejection and corre- 
sponds to enforcing the constraints Pg(y|x, z) = and 
P F (y, y|x, z) = 1— P A (y, z). For the reduced-flipping method, 
we wish to remove these constraints and minimize P F while 
imposing: 

P(x)S(x, y)P A (x, y)S(y, z)P F (y, y|x, z) 
+P(i)S(z, y)P A (z, y)S(y, x)P s (y|z, x) 
= P(x)5(x, y)P A (x, y)S(y, z) (1 - P A (y, z)) (8) 

noting that 

P F (y,y|x,z) 6 [0, 1 - P A (y, z)] (9) 

P s (y|z,x) G [0,l-P A (y,x)] (10) 

Thus, equation [8] considers the proposed forward transition 
y — > z given x and the proposed reverse transition y — > x 
given z. By considering the overall rejection rates, we can en- 
sure that the number of transitions to state y by the reduced- 
flipping algorithm (given by the left hand side of equation 

I 



P F (y,y|x,z) = 1 -P A (y,z) -P s (y|x,z) 



Equations [12] and [T3l which define the probabilities of 
self- and flipping-transitions for the reduced-flipping GHMC 
method, arc the central result of this manuscript. 

Once again, we wish to prove satisfaction of the balance 
condition by equating P* (y) , the probability of observing y 

I 



[5]) remains equivalent to the overall rejection rate of the pro- 
posed forward move (the right hand side of equation [8]) while 
minimizing the number of momentum flips resulting from each 
proposed transition. This method breaks the detailed balance 
condition for self-transitions, corresponding to the modified 
condition: 

P(x)T(x, y) = P(y)T(y, x), x / y (11) 

This reduced-flipping method applies only to the case in 
which the previous transition resulted from an accepted move. 
If the previous transition resulted from a rejected move, we 
will use the standard GHMC criteria (automatic momentum 
flipping). This limitation is discussed in more detail below. 

The probabilities given the above constraints are easily de- 
fined: 



(12) 
(13) 

I 

within our simulation, to the stationary distribution P(y). 
P*(y) receives contributions from all accepted transitions to 
y and from rejected candidates from the initial states y and 

y- 



min(l-P A (y,z),||f|£^f||(l-P A (y,x))) if x ^ y was an 
-Ps(y|x, z) — < accepted transition 



else 



P*(y) = / dxP(x)S(x,y)P A (x,y)+ f dx f dz [P(x)5(x, y)P A (x, y)5(y, z)P s (y|x, z) 
Jn Jn Jn 

+P(x)S(x, y)P A (x, y)5(y, z)P F (y, y|x, z)] 

= f dxP(y)5(y,x)P A (y,x)+ f dx f dz [P(x)S(x, y)P A (x, y)5(y, z)P s (y|x, z) 
Jn Jn Jn 

+P(i)5(z, y)P A (z, y)5(y, x)P F (y, y|z, x)] 

= P(y)r(y)+ I dx / dzP(x)5(x,y)P A (x,y)5(y,z)(l-P A (y,z)) 
Jn Jn 

= P(fMy)+ [ dx [ dzP(y)5(y,x)P A (y,x)S(y,z)(l-P A (y,z)) 
Jn Jn 

= P(?)r(y) + P(y) f dz5(y,z)(l - P A (y,z)) 
Jn 

= P(y)r(y) + P(y)(l-r(y))=P(y) (14) 

I 

As is demonstrated by equation [T3l the reduced-flipping balance condition given in equation [Til which cannot be ap- 
method reverts to the automatic-flipping algorithm if the pre- plied to transitions resulting from rejected moves. It is possi- 
vious transition was a rejected move. This is because the de- ble to remove this limitation by establishing self- and flipping- 
crease in momentum-flipping transitions relies on the detailed transition probabilities that depend on the path of rejections 
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Figure 2: C(t), the radius of gyration autocorrelation function, 
calculated using the VVVR integrator at both high ( 7 = 91 ps -1 ) 
and low ( 7 = 0.01 ps" 1 ) viscosity. All simulations were per- 
formed using a 10 fs time step and resulted in at least a 97% 
acceptance rate. Solid curves were obtained with no Metropo- 
lis test while dashed curves were obtained using the traditional 
GHMC automatic-flipping algorithm. As can be seen, use of 
this automatic-flipping algorithm significantly extends the time re- 
quired for system decorrelation for simulations performed with a 

small friction coefficient, while this effect is very small for Simula- zj. R,6SU.ltS 
tions performed with a high friction constant. 



Figure 3: C(t), the radius of gyration autocorrelation function, 
using a simulation performed with no Metropolis test as a stan- 
dard comparison. Simulations were performed using the VVVR 
integrator with 7 = 0.01 ps -1 . Dashed lines correspond to results 
obtained using the traditional, automatic-flipping algorithm while 
solid lines correspond to those obtained with the reduced-flipping 
algorithm. 



(and, thus, every rejected candidate state) taken following the 
last accepted transition. This modification adds significant 
complexity to the algorithm (the calculation of such proba- 
bilities scale as 0(N\) for N serial rejections) and is unlikely 
to add significant benefit to simulation mixing times. 

3 Methods 

All simulations were performed using a modified version of 
GROMACS-3.3.1 El and the MARTINI 2.1 forcefield 



[171 1 1 81 ] . Results were collected on an eight-residue repeat of 
polyglutaminc immersed in water at 400K using either the 
VVVR integrator II Oil or verlet integration with stochastic 
velocity rescaling [12j. Simulations using VVVR integration 
were performed over a range of values for the friction coef- 
ficient, 7 s {0.01, 1.0,91} ps -1 . Simulations were performed 
using nonbonded interactions that are switched at 0.9 nm and 
cut off at 1.2 nm in conjunction with a dispersive tail correc- 
tion [H]. 

As an indication of mixing times for the algorithms pre- 
sented in this work, we study the autocorrelation function for 
the radius of gyration, R g , for the polyglutaminc octamcr: 



C{5t) 



(R g (t)R g (t + 8t)) 



(15) 



We also analyze the integrated autocorrelation function, 
An = L C (6t) ddt. Samples were collected every 10 ps and 
statistical uncertainties for integrated autocorrelation func- 
tions were estimated using bootstrapping analysis. 



Figure [2] compares the autocorrelation functions for a poly- 
mer radius of gyration calculated from simulations using the 
traditional GHMC (automatic-flipping) algorithm and from 
those with no Metropolis test. These simulations were per- 
formed with a 10 fs time step, resulting in an acceptance rate 
of at least 97% for all GHMC simulations. As can be seen, 
the effect of the automatic-flipping algorithm is insignificant 
for simulations performed with the VVVR integrator and a 
friction constant of 7 = 91 ps -1 . The algorithm is shown 
to significantly extend the time required for system decor- 
relation, however, as the the friction constant is decreased. 
This effect is most pronounced for results obtained using the 
stochastic velocity rescaling thermostat. This is not surpris- 
ing, given that this integration scheme can be derived start- 
ing with Langevin integration and minimizing the thermo- 
stat's disturbance on the Hamiltonian trajectory [2(J. Table 
[T] gives the integrated autocorrelation function for these re- 
sults. Again, for the most pronounced case using the stochas- 
tic velocity rescaling thermostat, the An g value increases over 
an order of magnitude, from 0.07 ± .01 ns for a simulation 
containing no Metropolis test to 1.62 ± 0.02 ns for tradtional 
GHMC. 

For simulations performed using the VVVR integrator with 
low viscosity (7 = 0.01 ps -1 ), figure [3] displays autocorrelation 
functions for a polymer radius of gyration calculated using 
the automatic-flipping and reduced-flipping algorithms. The 
effect of the GHMC method on system dynamics is substan- 
tially mitigated by use of the reduced-flipping algorithm. For 
simulations performed with the stochastic velocity rescaling 
thermostat, An g — 0.22 ± 0.02 ns for the reduced-flipping 
algorithm. This is very close to the value obtained with no 
Metropolis test (An = 0.07 ± .01 ns), in contrast to the in- 
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Table 1: Comparison of the integrated autocorrelation function for simulations performed with no Metropolis test, traditional GHMC 
with automatic momentum flipping, and reduced-flipping GHMC. 



Timestep (fs) 

Algorithm 


10 


20 


30 


VVVR, 7 = 91 ps" 1 


Ar „(ns), no Metropolis test 
Ar „(ns), automatic-flipping 
Ar b (ns) , reduced-flipping 


4.61 ± 0.08 
5.41 ± 0.12 
4.99 ± 0.08 


6.10 ± 0.05 
5.50 ± 0.04 


10.23 ± 0.38 
10.08 ± 0.51 


VVVR, 7 = 1.0 ps- 1 


Ar (ns), no Metropolis test 
Ar (ns), automatic-flipping 
Ar 9 (ns) , reduced-flipping 


0.14 ± 0.01 
0.52 ± 0.01 
0.22 ± 0.01 


1.52 ± 0.06 
1.04 ±0.05 


11.50 ± 0.30 
7.11 ± 0.17 


VVVR, 7 = 0.01 ps" 1 


Ar (ns), no Metropolis test 
Ar^ (ns), automatic-flipping 
Ar 9 (ns) , reduced-flipping 


0.12 ± 0.01 
1.04 ± 0.01 
0.16 ± 0.01 


3.46 ± 0.09 
2.00 ± 0.04 


23.04 ± 0.31 
17.96 ± 0.42 


Stochastic velocity rescaling 


Ar (ns), no Metropolis test 
Ar (ns), automatic-flipping 
A Rg (ns) , reduced-flipping 


0.07 ± 0.01 
1.62 ± 0.02 
0.22 ± 0.02 


5.21 ± 0.09 
2.07 ± 0.04 


28.11 ± 0.19 
20.93 ± 0.16 



crease by an order of magnitude observed for the automatic- 
flipping algorithm (An g = 1.62 ± 0.02 ns). 

The benefit of reduced-flipping GHMC becomes negligible 
at large time steps and small acceptance rates (for dt = 30 fs, 
these GHMC simulations resulted in acceptance rates « 40%). 
Unsurprisingly, the reduced-flipping algorithm also demon- 
strates no major benefit for simulations performed at high 
viscosity (see table Q] for a comparison of automatic-flipping 
and reduced-flipping algorithms over all integration parame- 
ters). 

5 Conclusions 

The introduction of a Metropolis criterion into a molecu- 
lar dynamics simulation traditionally requires an undesirable 
transformation of the system momenta, either via complete 
randomization (HMC) or momentum flipping upon move re- 
jection (GHMC). This manuscript presents results obtained 
using a coarse-grained forccficld simulation of a polyglutamine 
octamer and demonstrates that concerns over the effect of 
these algorithms on simulation mixing times may be well jus- 
tified, depending on the particular method of integration. The 
two main conclusions of this paper are that (1) for certain 
sets of conditions, transformations of system momenta asso- 
ciated with HMC / GHMC algorithms significantly affect sim- 
ulation mixing times, and (2) our reduced-flipping algorithm 
presents a simple alternative to the traditional GHMC algo- 
rithm that significantly reduces this effect. These results have 
been demonstrated for both the VVVR integrator, equivalent 
to Langevin dynamics flfl ] , and the stochastic velocity rescal- 
ing thermostat, which can be derived as a form Langevin in- 
tegration by minimizing the thermostat's disturbance on the 
Hamiltonian trajectory [iol ]. 

The differences between results obtained with no Metropo- 
lis test, with traditional automatic- flipping GHMC, and with 
the reduced-flipping GHMC arc insignificant for simulations 
performed with a high friction coefficient but become quite 



pronounced for integration methods corresponding to low fric- 
tion or velocity collision rate. Surprisingly, these effects are 
significant even at small time steps associated with high ac- 
ceptance rates. The reduced-flipping algorithm presented in 
this work, a simple modification to the traditional GHMC 
method that rigorously satisfies the balance condition, is ca- 
pable of significantly reducing the effect that these automatic- 
flipping algorithms have on molecular kinetics and simulation 
mixing times. 
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7 Appendix 

A The selection probability for in- 
tegration steps using a stochastic 
velocity rescaling thermostat 



Thermostatting via stochastic velocity rescaling [12j is similar 
to the Bercndsen thermostat 2l| with an additional stochas- 
tic term that ensures proper canonical sampling. The posi- 
tions and velocities are propagated using standard integration 
methods. At iteration i, the kinetic energy (E^) is updated: 



K* = K 



(Ef-El)^ 



NfT 



(16) 



where E™ f 



^ is the average kinetic energy given the num- 



ber of degrees of freedom (Nf) and the inverse thermal energy 
(/3), the parameter t defines the relaxation timescale, and £ 



5 



is a Gaussian random variable of mean and standard devia- 
tion 1. To achieve this change in kinetic energy, all velocities 

are scaled by the factor 

We can calculate the Gaussian random variable £ necessary 
for the reverse transition: 



El 



>N 



2jE»EfAtT 



(17) 



Note that we must include the Jacobian factor accounting 
for the coordinate transformation from £ (a single random 
variable applied to the total kinetic energy) to the Cartesian 
space of the Nf momentum degrees of freedom. The ratio of 
selection probabilities for a transition (x — > y) is: 



S(x^y) 
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